options(width=108)
out=F
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(ggplot2)
library(plyr)
library(reshape)
library(matrixStats)
library(gridExtra)
library(lme4)
library(lmerTest)
library(BayesFactor)
#library(splines)
db <- '/home/egor/Dropbox/' # on Linux
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))
# theme for plotting:
alpha <- .7
w <- .56 # proportion width in group plots
themefy <- function(p) {
p <- p + theme_bw() +
theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
axis.text=element_text(size=8), axis.title=element_text(size=9),
legend.text=element_text(size=8), legend.title=element_text(size=9),
legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
legend.background = element_rect(fill='transparent'),
plot.title=element_text(face='bold'))
}
cc <- c('#e31a1c','#fdbf6f','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99')
xLab <- 'Target Peak Time (s)'
yLab <- 'Log Contrast Threshold'
# colLab <- expression(paste('\nTarget\nVelocity (', degree, '/s)', sep=''))
colLab <- expression(paste('Target Eccentricity (', degree, ')', sep=''))
# colLabType <- 'Mask Type'
yLim <- c(-1.5,-0.5)
dodge <- position_dodge(width=0.1)
sumFn <- function(ss, subjStr='subj', xStr='targTpeak', grpStr='targEcc'){
sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) #, .drop=F)
sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) # total mean across conditions per subj
sumSubj <- merge(sumSubj, sumSubjMn)
sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
sumSubj$seNorm <- NA
# sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
sumGrp$subj <- 'average'
sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
sumComb$se[is.na(sumComb$se)] <- 0
sumComb
}
plotAve <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
labs(x=xlab, y=ylab, colour=collab) + #ylim(0,1) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
plotIndiv <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
facet_wrap( ~ subj, ncol=4) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
# scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels=c('0','','0.5','','1'), limits=c(0,1)) +
labs(x=xlab, y=ylab, colour=collab) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
allDataDir <- paste(db,'Projects/mc/data_bv3/mcBv3_xvv',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('xvv',dataDirs)]
colsOfInt <- c('participant', 'dom', 'session', 'mcBv', 'targTpeak', 'targXoff2',
'targV', 'stairStart', 'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
print(curDir)
curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
curDf <- curDf[,colsOfInt]
curDf$cond <- substr(curDir,19,22)
subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
subjStairs <- subjStairs[grep('.tsv',subjStairs)]
for(curStairFN in subjStairs){
curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
curRevs <- read.table(curStair, skip=1, nrows=1)
curIntn <- read.table(curStair, skip=4, nrows=2)
curInfo <- readLines(curStair)
curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1],
sess=curDf$session[1], stairStart=curInfo[33],
mcBv=curInfo[41], targTpeak=curInfo[43],
targEcc=curInfo[37], targV=curInfo[31])
curDfRevs <- cbind(curInfoCols, curRevs[,2:9])
nTrials <- ncol(curIntn)-1
curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
curDfIntn$trial <- seq(1,nTrials)
curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
rownames(curDfIntn) <- NULL
dfRevs <- rbind(dfRevs, curDfRevs)
dfIntn <- rbind(dfIntn, curDfIntn)
}
df <- rbind(df, curDf)
}
## [1] "mc2_tgT-mcBv3_xvv-cent_p1_s1_2017-10-05_1028"
## [1] "mc2_tgT-mcBv3_xvv-cent_p12_s1_2017-10-10_1604"
## [1] "mc2_tgT-mcBv3_xvv-cent_p13_s1_2017-10-10_1653"
## [1] "mc2_tgT-mcBv3_xvv-cent_p15_s1_2017-10-12_1330"
## [1] "mc2_tgT-mcBv3_xvv-cent_p5_s2_2017-10-06_1627"
## [1] "mc2_tgT-mcBv3_xvv-cent_p7_s2_2017-10-06_1432"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p10_s1_2017-10-10_1033"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p11_s2_2017-10-10_1425"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p14_s2_2017-10-13_1137"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p16_s1_2017-10-13_1615"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p2_s1_2017-10-05_1234"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p6_s1_2017-10-06_1040"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p9_s1_2017-10-06_1549"
## [1] "mc2_tgT-mcBv3_xvv-peri_p13_s2_2017-10-12_1608"
## [1] "mc2_tgT-mcBv3_xvv-peri_p15_s2_2017-10-12_1410"
## [1] "mc2_tgT-mcBv3_xvv-peri_p4_s1_2017-10-05_1539"
## [1] "mc2_tgT-mcBv3_xvv-peri_p5_s1_2017-10-05_1640"
## [1] "mc2_tgT-mcBv3_xvv-peri_p7_s1_2017-10-06_1245"
## [1] "mc2_tgT-mcBv3_xvv-stat_p11_s1_2017-10-10_1247"
## [1] "mc2_tgT-mcBv3_xvv-stat_p14_s1_2017-10-12_1240"
## [1] "mc2_tgT-mcBv3_xvv-stat_p2_s2_2017-10-05_1425"
## [1] "mc2_tgT-mcBv3_xvv-stat_p3_s1_2017-10-05_1342"
## [1] "mc2_tgT-mcBv3_xvv-stat_p6_s2_2017-10-12_1035"
## [1] "mc2_tgT-mcBv3_xvv-stat_p8_s1_2017-10-06_1357"
## [1] "mc2_tgT-mcBv3_xvv-stat_p9_s2_2017-10-10_1348"
ds <- rename(df, c(participant='subj', session='sess', meanRev6='thresh',
mcBv='maskV', targXoff2='targEcc'))
ds$targEcc <- round(ds$targEcc / 35,1)
ds$targV <- round(ds$targV / 3.5,1)
ds$maskType <- ''
ds$maskType[ds$maskV==0.01] <- 'static'
ds$maskType[ds$maskV==0.6] <- 'slow'
ds$maskType[ds$maskV==9.6] <- 'fast'
ds$maskV <- round(ds$maskV * 60 / 35, 1)
ds$maskV[ds$maskV<0.05] <- 0
dfIntn$maskV <- round(as.numeric(levels(dfIntn$mcBv))[dfIntn$mcBv] * 60 / 35, 1)
dfIntn$maskV[dfIntn$maskV<0.05] <- 0
ds$condSplit <- ''
ds$condSplit[ds$cond=='stat' | ds$cond=='dyna'] <- 'v'
ds$condSplit[ds$cond=='cent' | ds$cond=='peri'] <- 'ecc'
head(ds)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit
## 1 1 1 1 0.0 0.5 0.8 0.0 -2 -0.9583333 cent static ecc
## 2 1 1 1 1.0 0.5 0.8 4.6 -2 -0.7166667 cent slow ecc
## 3 1 1 1 16.5 1.0 0.8 0.0 -2 -1.6083333 cent fast ecc
## 4 1 1 1 1.0 1.5 0.8 4.6 0 -0.6583333 cent slow ecc
## 5 1 1 1 1.0 1.0 0.8 0.0 -2 -0.8750000 cent slow ecc
## 6 1 1 1 16.5 0.5 0.8 0.0 -2 -1.5916667 cent fast ecc
thresh <- ddply(ds, .(subj,dom,sess,maskV,maskType,targTpeak,targEcc,targV,
cond,condSplit), summarise, threshMean = mean(thresh))
head(thresh)
## subj dom sess maskV maskType targTpeak targEcc targV cond condSplit threshMean
## 1 1 1 1 0 static 0.5 0.8 0.0 cent ecc -0.8916667
## 2 1 1 1 0 static 0.5 0.8 4.6 cent ecc -1.2208333
## 3 1 1 1 0 static 1.0 0.8 0.0 cent ecc -1.1041667
## 4 1 1 1 0 static 1.0 0.8 4.6 cent ecc -1.2333333
## 5 1 1 1 0 static 1.5 0.8 0.0 cent ecc -1.0875000
## 6 1 1 1 0 static 1.5 0.8 4.6 cent ecc -1.2375000
ss <- sumFn(thresh[(thresh$maskType=='static' & thresh$targV==0),])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==0 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='static' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==0 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='slow' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==1 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='slow' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==1 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==16.5 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==16.5 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
grid.arrange(pTargStatMaskStat + #theme(legend.position='none') +
ggtitle('a: static mask, static target'),
pTargDynaMaskStat + ggtitle('d: static mask, dynamic target'),
pTargStatMaskSlow + ggtitle('b: slow mask, static target'),
pTargDynaMaskSlow + ggtitle('e: slow mask, dynamic target'),
pTargStatMaskFast + ggtitle('c: fast mask, static target'),
pTargDynaMaskFast + ggtitle('f: fast mask, dynamic target'),
ncol=2, nrow=3) #widths=c((1-w)/2,w/2))
cent <- function(v){
v <- apply(v,2,function(x){
x <- x - mean(unique(x),na.rm=T)
x <- x / max(x)
})
}
dsc <- ds
centCols <- c('dom','targTpeak','stairStart','targEcc','maskV','targV','sess')
dsc[,centCols] <- cent(dsc[,centCols])
# dsc$maskV <- (dsc$maskV / max(dsc$maskV)) #*2 - 1
dsc$maskV_c <- (dsc$maskV + 1) / 2
dsc$targV_c <- (dsc$targV + 1) / 2
dsc$targEcc_c <- round((dsc$targEcc + 1) / 2,0)
dsc$targTpeak_c <- (dsc$targTpeak + 1) / 2
# dsc$subj <- as.factor(dsc$subj)
# dsc$maskType <- ordered(ds$maskType, levels=c('static','slow','fast'))
dsc$maskType <- factor(dsc$maskType, c('static','slow','fast'))
head(dsc)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit maskV_c
## 1 1 1 -1 -0.546875 -1 -1 -1 -1 -0.9583333 cent static ecc 0.2265625
## 2 1 1 -1 -0.453125 -1 -1 1 -1 -0.7166667 cent slow ecc 0.2734375
## 3 1 1 -1 1.000000 0 -1 -1 -1 -1.6083333 cent fast ecc 1.0000000
## 4 1 1 -1 -0.453125 1 -1 1 1 -0.6583333 cent slow ecc 0.2734375
## 5 1 1 -1 -0.453125 0 -1 -1 -1 -0.8750000 cent slow ecc 0.2734375
## 6 1 1 -1 1.000000 -1 -1 -1 -1 -1.5916667 cent fast ecc 1.0000000
## targV_c targEcc_c targTpeak_c
## 1 0 0 0.0
## 2 1 0 0.0
## 3 0 0 0.5
## 4 1 0 1.0
## 5 0 0 0.5
## 6 0 0 0.0
# ds$targTpeak_c <- ds$targTpeak - 0.5 # [.5, 1, 1.5] -> [0, .5, 1]
# ds$targEcc_c <- cent(ds$targEcc)
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c +
(1|subj), data=dsc[dsc$targV==-1,]))
## estm se df tVal pVal star
## (Intercept) -0.993 0.064 23.951 -15.595 0.000 ***
## dom -0.020 0.053 11.970 -0.376 0.714
## stairStart 0.029 0.009 423.116 3.124 0.002 **
## sess 0.057 0.016 420.799 3.568 0.000 ***
## maskTypeslow 0.173 0.050 423.116 3.447 0.001 ***
## maskTypefast -0.330 0.050 423.116 -6.575 0.000 ***
## targTpeak_c -0.077 0.055 423.116 -1.404 0.161
## targEcc_c 0.083 0.052 423.964 1.610 0.108
## maskTypeslow:targTpeak_c 0.028 0.078 423.116 0.362 0.717
## maskTypefast:targTpeak_c 0.019 0.078 423.116 0.243 0.808
## maskTypeslow:targEcc_c -0.293 0.073 423.116 -4.045 0.000 ***
## maskTypefast:targEcc_c -0.061 0.073 423.116 -0.838 0.403
## targTpeak_c:targEcc_c 0.037 0.079 423.116 0.465 0.642
## maskTypeslow:targTpeak_c:targEcc_c 0.035 0.112 423.116 0.308 0.758
## maskTypefast:targTpeak_c:targEcc_c -0.003 0.112 423.116 -0.023 0.982
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c *
targV_c + (1|subj), data=dsc))
## estm se df tVal pVal star
## (Intercept) -1.027 0.067 24.978 -15.241 0.000 ***
## dom 0.030 0.058 13.926 0.507 0.620
## stairStart 0.029 0.006 858.941 4.593 0.000 ***
## sess -0.004 0.008 865.348 -0.545 0.586
## maskTypeslow 0.173 0.049 858.941 3.565 0.000 ***
## maskTypefast -0.330 0.049 858.941 -6.798 0.000 ***
## targTpeak_c -0.077 0.053 858.941 -1.451 0.147
## targEcc_c 0.085 0.050 859.259 1.698 0.090 .
## targV_c 0.036 0.049 859.390 0.735 0.463
## maskTypeslow:targTpeak_c 0.028 0.075 858.941 0.375 0.708
## maskTypefast:targTpeak_c 0.019 0.075 858.941 0.251 0.802
## maskTypeslow:targEcc_c -0.293 0.070 858.941 -4.182 0.000 ***
## maskTypefast:targEcc_c -0.061 0.070 858.941 -0.866 0.387
## targTpeak_c:targEcc_c 0.037 0.077 858.941 0.481 0.630
## maskTypeslow:targV_c 0.064 0.069 858.941 0.938 0.349
## maskTypefast:targV_c 0.490 0.069 858.941 7.138 0.000 ***
## targTpeak_c:targV_c 0.038 0.075 858.941 0.502 0.615
## targEcc_c:targV_c 0.057 0.070 858.942 0.817 0.414
## maskTypeslow:targTpeak_c:targEcc_c 0.035 0.109 858.941 0.319 0.750
## maskTypefast:targTpeak_c:targEcc_c -0.003 0.109 858.941 -0.024 0.981
## maskTypeslow:targTpeak_c:targV_c -0.014 0.106 858.941 -0.132 0.895
## maskTypefast:targTpeak_c:targV_c 0.032 0.106 858.941 0.304 0.761
## maskTypeslow:targEcc_c:targV_c 0.134 0.099 858.941 1.352 0.177
## maskTypefast:targEcc_c:targV_c -0.030 0.099 858.941 -0.299 0.765
## targTpeak_c:targEcc_c:targV_c -0.052 0.109 858.941 -0.476 0.634
## maskTypeslow:targTpeak_c:targEcc_c:targV_c -0.016 0.154 858.941 -0.105 0.917
## maskTypefast:targTpeak_c:targEcc_c:targV_c -0.010 0.154 858.941 -0.066 0.947